#include "kernel.h"


float kernel(glm::vec3 t, float dist) {
    float d;

    d = dist * dist - glm::dot(t, t);

    if (dist < glm::length(t))
        return 0.f;

    /* (315 * (dist^2 - glm::dot(t, t))^3) / (64 * pi * dist^9) */
    return 315.f * (d * d * d) / (64.f * PI * 
        dist * dist * dist * 
        dist * dist * dist *
        dist * dist * dist);
}

glm::vec3 grad_kernel(glm::vec3 t, float dist) {
    float d;

    d = glm::length(t);
    d = dist * dist - d * d; 

    if (dist < glm::length(t))
        return 0.f * t;

    return t * 
        (-945.f * (d * d) / (32.f * PI * 
        dist * dist * dist * 
        dist * dist * dist *
        dist * dist * dist));
}

float laplacian_kernel(glm::vec3 t, float dist) {
    float d;

    d = glm::length(t);

    if (d > dist)
        return 0.f;

    /* (315 * (dist^2 - glm::dot(t, t))^3) / (64 * pi * dist^9) */
    return (-945.f * (dist * dist - d * d) * 
     (3 * dist * dist - 7 * d * d) / 
     (32.f * PI * 
     dist * dist * dist * 
     dist * dist * dist *
     dist * dist * dist));
}

/* pressure kernels */
float pressure_kernel(glm::vec3 t, float dist) {
    float d;

    d = dist - glm::length(t);

    if (d < 0.f)
        return 0.f;

    /* (15 * (dist - glm::length(t))^3) / (pi * dist^6) */
    return 15.f * (d * d * d) / (PI * 
        dist * dist * dist *
        dist * dist * dist);
}

glm::vec3 limit_pressure(glm::vec3 t, float dist) {
    glm::vec3 ret;
    float dist_6 = dist * dist * dist * dist * dist * dist;

    if (
     t.x > 0.001 || t.x < -0.001 ||
     t.y > 0.001 || t.y < -0.001 ||
     t.z > 0.001 || t.z < -0.001)
        printf("got here\n");
    ret.x = 
     t.x < 0.001 && t.x > 0 ? -45.f / PI / dist_6 :
     t.x > -0.001 && t.x <= 0 ? 45.f / PI / dist_6 :
     45.f * t.x * dist / PI / dist_6;
    ret.y = 
     t.y < 0.001 && t.y > 0 ? -45.f / PI / dist_6 :
     t.y > -0.001 && t.y <= 0 ? 45.f / PI / dist_6 :
     45.f * t.x * dist / PI / dist_6;
    ret.z = 
     t.z < 0.001 && t.z > 0 ? -45.f / PI / dist_6 :
     t.z > -0.001 && t.z <= 0 ? 45.f / PI / dist_6 :
     45.f * t.x * dist / PI / dist_6;
    
    return ret;
}

glm::vec3 grad_pressure(glm::vec3 t, float dist) {
    float d;
    float t_dot_t;

    d = dist - glm::length(t);

    /*if (d < 0.f)
        return 0.f * t;
    else if (glm::length(t) <= 0.001) {
        return limit_pressure(t, dist);
    }*/

    return glm::normalize(t) * (-45.f * (d * d) / (PI * 
        dist * dist * dist *
        dist * dist * dist));
}

float pressure_laplacian(glm::vec3 t, float dist) {
    float d;

    d = glm::length(t);

    /*if (d > dist)
        return 0.f;

    if (d < 0.0001)
        return -10000.f;*/

    return (-90.f * (dist - d) * (dist - 2.f * d) / 
     (PI * d *
     dist * dist * dist * 
     dist * dist * dist *
     dist * dist * dist));
}

/* viscosity kernels */
float viscosity_kernel(glm::vec3 t, float dist) {
    float d;

    d = glm::length(t);

    if (d > dist)
        return 0.f;

    if (d < 0.0001)
        return 1000.f;

    return 15.f * 
     (-(d * d * d / (2.f * dist * dist * dist)) + 
     d * d / dist / dist + dist / (2.f * d) - 1) / 
     (2.f * PI * dist * dist * dist);
}

glm::vec3 viscosity_lim(glm::vec3 t, float dist) {
    glm::vec3 ret;

    if (t.x <= 0.f && t.x > -0.07) ret.x = 100.f;
    else if (t.x > 0.f && t.x < 0.07) ret.x = -100.f;
    else ret.x = 0.0f;

    if (t.y <= 0.f && t.y > -0.07) ret.y = 100.f;
    else if (t.y > 0.f && t.y < 0.07) ret.y = -100.f;
    else ret.y = 0.0f;

    if (t.z <= 0.f && t.z > -0.07) ret.z = 100.f;
    else if (t.z > 0.f && t.z < 0.07) ret.z = -100.f;
    else ret.z = 0.0f;

    return ret;
}

glm::vec3 grad_viscosity(glm::vec3 t, float dist) {
    float d;

    d = glm::length(t);

    /*if (d > dist)
        return glm::vec3(0.f);

    if (d < 0.0001) {
        return viscosity_lim(t, dist);
    }*/

    return 15.f * t *
     (-(d * 3.f / (2.f * dist * dist * dist)) + 
     2.f / dist / dist + dist / (2.f * d * d * d)) / 
     (2.f * PI * dist * dist * dist);
}

float laplacian_viscosity(glm::vec3 t, float dist) {
    float d;

    d = dist - glm::length(t);

    if (d < 0.f)
        return 0.f;
    
    return (45.f * (d) / (PI * 
        dist * dist * dist *
        dist * dist * dist));
}
